TCT/TNF joint analysis to identify MKi response markers

Index

• Gene level analysis
• Transcript level analysis

Library loading and set up

suppressMessages(
  c(library(DESeq2),
    library(tximport),
    library(gridExtra),
    library(ensembldb),
    library(EnsDb.Mmusculus.v79),
    library(grid),
    library(ggplot2),
    library(lattice),
    library(reshape),
    library(mixOmics),
    library(gplots),
    library(RColorBrewer),
    library(readr),
    library(dplyr),
    library(VennDiagram),
    library(org.Mm.eg.db), 
    library(pathview),
    library(AnnotationDbi),
    library(clusterProfiler),
    library(fgsea))
)

##   [1] "DESeq2"               "SummarizedExperiment" "DelayedArray"        
##   [4] "matrixStats"          "Biobase"              "GenomicRanges"       
##   [7] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
##  [10] "BiocGenerics"         "parallel"             "stats4"              
##  [13] "stats"                "graphics"             "grDevices"           
##  [16] "utils"                "datasets"             "methods"             
##  [19] "base"                 "tximport"             "DESeq2"              
##  [22] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
##  [25] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
##  [28] "IRanges"              "S4Vectors"            "BiocGenerics"        
##  [31] "parallel"             "stats4"               "stats"               
##  [34] "graphics"             "grDevices"            "utils"               
##  [37] "datasets"             "methods"              "base"                
##  [40] "gridExtra"            "tximport"             "DESeq2"              
##  [43] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
##  [46] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
##  [49] "IRanges"              "S4Vectors"            "BiocGenerics"        
##  [52] "parallel"             "stats4"               "stats"               
##  [55] "graphics"             "grDevices"            "utils"               
##  [58] "datasets"             "methods"              "base"                
##  [61] "ensembldb"            "AnnotationFilter"     "GenomicFeatures"     
##  [64] "AnnotationDbi"        "gridExtra"            "tximport"            
##  [67] "DESeq2"               "SummarizedExperiment" "DelayedArray"        
##  [70] "matrixStats"          "Biobase"              "GenomicRanges"       
##  [73] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
##  [76] "BiocGenerics"         "parallel"             "stats4"              
##  [79] "stats"                "graphics"             "grDevices"           
##  [82] "utils"                "datasets"             "methods"             
##  [85] "base"                 "EnsDb.Mmusculus.v79"  "ensembldb"           
##  [88] "AnnotationFilter"     "GenomicFeatures"      "AnnotationDbi"       
##  [91] "gridExtra"            "tximport"             "DESeq2"              
##  [94] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
##  [97] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
## [100] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [103] "parallel"             "stats4"               "stats"               
## [106] "graphics"             "grDevices"            "utils"               
## [109] "datasets"             "methods"              "base"                
## [112] "grid"                 "EnsDb.Mmusculus.v79"  "ensembldb"           
## [115] "AnnotationFilter"     "GenomicFeatures"      "AnnotationDbi"       
## [118] "gridExtra"            "tximport"             "DESeq2"              
## [121] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
## [124] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
## [127] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [130] "parallel"             "stats4"               "stats"               
## [133] "graphics"             "grDevices"            "utils"               
## [136] "datasets"             "methods"              "base"                
## [139] "ggplot2"              "grid"                 "EnsDb.Mmusculus.v79" 
## [142] "ensembldb"            "AnnotationFilter"     "GenomicFeatures"     
## [145] "AnnotationDbi"        "gridExtra"            "tximport"            
## [148] "DESeq2"               "SummarizedExperiment" "DelayedArray"        
## [151] "matrixStats"          "Biobase"              "GenomicRanges"       
## [154] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
## [157] "BiocGenerics"         "parallel"             "stats4"              
## [160] "stats"                "graphics"             "grDevices"           
## [163] "utils"                "datasets"             "methods"             
## [166] "base"                 "lattice"              "ggplot2"             
## [169] "grid"                 "EnsDb.Mmusculus.v79"  "ensembldb"           
## [172] "AnnotationFilter"     "GenomicFeatures"      "AnnotationDbi"       
## [175] "gridExtra"            "tximport"             "DESeq2"              
## [178] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
## [181] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
## [184] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [187] "parallel"             "stats4"               "stats"               
## [190] "graphics"             "grDevices"            "utils"               
## [193] "datasets"             "methods"              "base"                
## [196] "reshape"              "lattice"              "ggplot2"             
## [199] "grid"                 "EnsDb.Mmusculus.v79"  "ensembldb"           
## [202] "AnnotationFilter"     "GenomicFeatures"      "AnnotationDbi"       
## [205] "gridExtra"            "tximport"             "DESeq2"              
## [208] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
## [211] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
## [214] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [217] "parallel"             "stats4"               "stats"               
## [220] "graphics"             "grDevices"            "utils"               
## [223] "datasets"             "methods"              "base"                
## [226] "mixOmics"             "MASS"                 "reshape"             
## [229] "lattice"              "ggplot2"              "grid"                
## [232] "EnsDb.Mmusculus.v79"  "ensembldb"            "AnnotationFilter"    
## [235] "GenomicFeatures"      "AnnotationDbi"        "gridExtra"           
## [238] "tximport"             "DESeq2"               "SummarizedExperiment"
## [241] "DelayedArray"         "matrixStats"          "Biobase"             
## [244] "GenomicRanges"        "GenomeInfoDb"         "IRanges"             
## [247] "S4Vectors"            "BiocGenerics"         "parallel"            
## [250] "stats4"               "stats"                "graphics"            
## [253] "grDevices"            "utils"                "datasets"            
## [256] "methods"              "base"                 "gplots"              
## [259] "mixOmics"             "MASS"                 "reshape"             
## [262] "lattice"              "ggplot2"              "grid"                
## [265] "EnsDb.Mmusculus.v79"  "ensembldb"            "AnnotationFilter"    
## [268] "GenomicFeatures"      "AnnotationDbi"        "gridExtra"           
## [271] "tximport"             "DESeq2"               "SummarizedExperiment"
## [274] "DelayedArray"         "matrixStats"          "Biobase"             
## [277] "GenomicRanges"        "GenomeInfoDb"         "IRanges"             
## [280] "S4Vectors"            "BiocGenerics"         "parallel"            
## [283] "stats4"               "stats"                "graphics"            
## [286] "grDevices"            "utils"                "datasets"            
## [289] "methods"              "base"                 "RColorBrewer"        
## [292] "gplots"               "mixOmics"             "MASS"                
## [295] "reshape"              "lattice"              "ggplot2"             
## [298] "grid"                 "EnsDb.Mmusculus.v79"  "ensembldb"           
## [301] "AnnotationFilter"     "GenomicFeatures"      "AnnotationDbi"       
## [304] "gridExtra"            "tximport"             "DESeq2"              
## [307] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
## [310] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
## [313] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [316] "parallel"             "stats4"               "stats"               
## [319] "graphics"             "grDevices"            "utils"               
## [322] "datasets"             "methods"              "base"                
## [325] "readr"                "RColorBrewer"         "gplots"              
## [328] "mixOmics"             "MASS"                 "reshape"             
## [331] "lattice"              "ggplot2"              "grid"                
## [334] "EnsDb.Mmusculus.v79"  "ensembldb"            "AnnotationFilter"    
## [337] "GenomicFeatures"      "AnnotationDbi"        "gridExtra"           
## [340] "tximport"             "DESeq2"               "SummarizedExperiment"
## [343] "DelayedArray"         "matrixStats"          "Biobase"             
## [346] "GenomicRanges"        "GenomeInfoDb"         "IRanges"             
## [349] "S4Vectors"            "BiocGenerics"         "parallel"            
## [352] "stats4"               "stats"                "graphics"            
## [355] "grDevices"            "utils"                "datasets"            
## [358] "methods"              "base"                 "dplyr"               
## [361] "readr"                "RColorBrewer"         "gplots"              
## [364] "mixOmics"             "MASS"                 "reshape"             
## [367] "lattice"              "ggplot2"              "grid"                
## [370] "EnsDb.Mmusculus.v79"  "ensembldb"            "AnnotationFilter"    
## [373] "GenomicFeatures"      "AnnotationDbi"        "gridExtra"           
## [376] "tximport"             "DESeq2"               "SummarizedExperiment"
## [379] "DelayedArray"         "matrixStats"          "Biobase"             
## [382] "GenomicRanges"        "GenomeInfoDb"         "IRanges"             
## [385] "S4Vectors"            "BiocGenerics"         "parallel"            
## [388] "stats4"               "stats"                "graphics"            
## [391] "grDevices"            "utils"                "datasets"            
## [394] "methods"              "base"                 "VennDiagram"         
## [397] "futile.logger"        "dplyr"                "readr"               
## [400] "RColorBrewer"         "gplots"               "mixOmics"            
## [403] "MASS"                 "reshape"              "lattice"             
## [406] "ggplot2"              "grid"                 "EnsDb.Mmusculus.v79" 
## [409] "ensembldb"            "AnnotationFilter"     "GenomicFeatures"     
## [412] "AnnotationDbi"        "gridExtra"            "tximport"            
## [415] "DESeq2"               "SummarizedExperiment" "DelayedArray"        
## [418] "matrixStats"          "Biobase"              "GenomicRanges"       
## [421] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
## [424] "BiocGenerics"         "parallel"             "stats4"              
## [427] "stats"                "graphics"             "grDevices"           
## [430] "utils"                "datasets"             "methods"             
## [433] "base"                 "org.Mm.eg.db"         "VennDiagram"         
## [436] "futile.logger"        "dplyr"                "readr"               
## [439] "RColorBrewer"         "gplots"               "mixOmics"            
## [442] "MASS"                 "reshape"              "lattice"             
## [445] "ggplot2"              "grid"                 "EnsDb.Mmusculus.v79" 
## [448] "ensembldb"            "AnnotationFilter"     "GenomicFeatures"     
## [451] "AnnotationDbi"        "gridExtra"            "tximport"            
## [454] "DESeq2"               "SummarizedExperiment" "DelayedArray"        
## [457] "matrixStats"          "Biobase"              "GenomicRanges"       
## [460] "GenomeInfoDb"         "IRanges"              "S4Vectors"           
## [463] "BiocGenerics"         "parallel"             "stats4"              
## [466] "stats"                "graphics"             "grDevices"           
## [469] "utils"                "datasets"             "methods"             
## [472] "base"                 "pathview"             "org.Hs.eg.db"        
## [475] "org.Mm.eg.db"         "VennDiagram"          "futile.logger"       
## [478] "dplyr"                "readr"                "RColorBrewer"        
## [481] "gplots"               "mixOmics"             "MASS"                
## [484] "reshape"              "lattice"              "ggplot2"             
## [487] "grid"                 "EnsDb.Mmusculus.v79"  "ensembldb"           
## [490] "AnnotationFilter"     "GenomicFeatures"      "AnnotationDbi"       
## [493] "gridExtra"            "tximport"             "DESeq2"              
## [496] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
## [499] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
## [502] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [505] "parallel"             "stats4"               "stats"               
## [508] "graphics"             "grDevices"            "utils"               
## [511] "datasets"             "methods"              "base"                
## [514] "pathview"             "org.Hs.eg.db"         "org.Mm.eg.db"        
## [517] "VennDiagram"          "futile.logger"        "dplyr"               
## [520] "readr"                "RColorBrewer"         "gplots"              
## [523] "mixOmics"             "MASS"                 "reshape"             
## [526] "lattice"              "ggplot2"              "grid"                
## [529] "EnsDb.Mmusculus.v79"  "ensembldb"            "AnnotationFilter"    
## [532] "GenomicFeatures"      "AnnotationDbi"        "gridExtra"           
## [535] "tximport"             "DESeq2"               "SummarizedExperiment"
## [538] "DelayedArray"         "matrixStats"          "Biobase"             
## [541] "GenomicRanges"        "GenomeInfoDb"         "IRanges"             
## [544] "S4Vectors"            "BiocGenerics"         "parallel"            
## [547] "stats4"               "stats"                "graphics"            
## [550] "grDevices"            "utils"                "datasets"            
## [553] "methods"              "base"                 "clusterProfiler"     
## [556] "pathview"             "org.Hs.eg.db"         "org.Mm.eg.db"        
## [559] "VennDiagram"          "futile.logger"        "dplyr"               
## [562] "readr"                "RColorBrewer"         "gplots"              
## [565] "mixOmics"             "MASS"                 "reshape"             
## [568] "lattice"              "ggplot2"              "grid"                
## [571] "EnsDb.Mmusculus.v79"  "ensembldb"            "AnnotationFilter"    
## [574] "GenomicFeatures"      "AnnotationDbi"        "gridExtra"           
## [577] "tximport"             "DESeq2"               "SummarizedExperiment"
## [580] "DelayedArray"         "matrixStats"          "Biobase"             
## [583] "GenomicRanges"        "GenomeInfoDb"         "IRanges"             
## [586] "S4Vectors"            "BiocGenerics"         "parallel"            
## [589] "stats4"               "stats"                "graphics"            
## [592] "grDevices"            "utils"                "datasets"            
## [595] "methods"              "base"                 "fgsea"               
## [598] "clusterProfiler"      "pathview"             "org.Hs.eg.db"        
## [601] "org.Mm.eg.db"         "VennDiagram"          "futile.logger"       
## [604] "dplyr"                "readr"                "RColorBrewer"        
## [607] "gplots"               "mixOmics"             "MASS"                
## [610] "reshape"              "lattice"              "ggplot2"             
## [613] "grid"                 "EnsDb.Mmusculus.v79"  "ensembldb"           
## [616] "AnnotationFilter"     "GenomicFeatures"      "AnnotationDbi"       
## [619] "gridExtra"            "tximport"             "DESeq2"              
## [622] "SummarizedExperiment" "DelayedArray"         "matrixStats"         
## [625] "Biobase"              "GenomicRanges"        "GenomeInfoDb"        
## [628] "IRanges"              "S4Vectors"            "BiocGenerics"        
## [631] "parallel"             "stats4"               "stats"               
## [634] "graphics"             "grDevices"            "utils"               
## [637] "datasets"             "methods"              "base"

Gene level analysis

Identify inflammation related DGE

First I need to identify DEG in inflammation in TCT and TNF models.

# DGE in inflammation using I_V and NI_V comparison
tct_inf <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Significant genes_I_V vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_6 = col_double(),
##   NI_V_7 = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_V_6 = col_double()
## )

tct_inf_DGE <- unique(tct_inf$X1)
tct_inf_DGE_up <- tct_inf$X1[tct_inf$log2FoldChange > 0]
tct_inf_DGE_down <- tct_inf$X1[tct_inf$log2FoldChange < 0]

tnf_inf <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Significant genes_I_V vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_3 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_5 = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double()
## )

tnf_inf_DGE <- unique(tnf_inf$X1)
tnf_inf_DGE_up <- tnf_inf$X1[tnf_inf$log2FoldChange > 0]
tnf_inf_DGE_down <- tnf_inf$X1[tnf_inf$log2FoldChange < 0]

# Find out how big is the overlap
## Total overlap
inf_overlap <-intersect(tct_inf_DGE, tnf_inf_DGE)

grid.newpage()
draw.pairwise.venn(length(tct_inf_DGE),
                   length(tnf_inf_DGE),
                   length(inf_overlap),
                   catergory <- c("TCT inflammation DGE",
                                  "TNF inflammation DGE"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.11], polygon[GRID.polygon.12], polygon[GRID.polygon.13], polygon[GRID.polygon.14], text[GRID.text.15], text[GRID.text.16], text[GRID.text.17], text[GRID.text.18], text[GRID.text.19])

pdf("TCT_TNF_inflammation_DGE_overlap.pdf",
           width = 14,
           height = 10)
draw.pairwise.venn(length(tct_inf_DGE),
                   length(tnf_inf_DGE),
                   length(inf_overlap),
                   catergory <- c("TCT inflammation DGE",
                                  "TNF inflammation DGE"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.20], polygon[GRID.polygon.21], polygon[GRID.polygon.22], polygon[GRID.polygon.23], text[GRID.text.24], text[GRID.text.25], text[GRID.text.26], text[GRID.text.27], text[GRID.text.28])

dev.off()

## quartz_off_screen 
##                 2

## Overlap in up-regulated genes
inf_overlap_up <-intersect(tct_inf_DGE_up, tnf_inf_DGE_up)

grid.newpage()
draw.pairwise.venn(length(tct_inf_DGE_up),
                   length(tnf_inf_DGE_up),
                   length(inf_overlap_up),
                   catergory <- c("TCT inflammation Up-DGE",
                                  "TNF inflammation Up-DGE"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.29], polygon[GRID.polygon.30], polygon[GRID.polygon.31], polygon[GRID.polygon.32], text[GRID.text.33], text[GRID.text.34], text[GRID.text.35], text[GRID.text.36], text[GRID.text.37])

pdf("TCT_TNF_inflammation_DGE_overlap_up.pdf",
           width = 14,
           height = 10)
draw.pairwise.venn(length(tct_inf_DGE_up),
                   length(tnf_inf_DGE_up),
                   length(inf_overlap_up),
                   catergory <- c("TCT inflammation Up-DGE",
                                  "TNF inflammation Up-DGE"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.38], polygon[GRID.polygon.39], polygon[GRID.polygon.40], polygon[GRID.polygon.41], text[GRID.text.42], text[GRID.text.43], text[GRID.text.44], text[GRID.text.45], text[GRID.text.46])

dev.off()

## quartz_off_screen 
##                 2

## Overlap in down-regulated genes
inf_overlap_down <-intersect(tct_inf_DGE_down, tnf_inf_DGE_down)

grid.newpage()
draw.pairwise.venn(length(tct_inf_DGE_down),
                   length(tnf_inf_DGE_down),
                   length(inf_overlap_down),
                   catergory <- c("TCT inflammation Down-DGE",
                                  "TNF inflammation Down-DGE"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.47], polygon[GRID.polygon.48], polygon[GRID.polygon.49], polygon[GRID.polygon.50], text[GRID.text.51], text[GRID.text.52], lines[GRID.lines.53], text[GRID.text.54], lines[GRID.lines.55], text[GRID.text.56], text[GRID.text.57])

pdf("TCT_TNF_inflammation_DGE_overlap_down.pdf",
           width = 14,
           height = 10)
draw.pairwise.venn(length(tct_inf_DGE_down),
                   length(tnf_inf_DGE_down),
                   length(inf_overlap_down),
                   catergory <- c("TCT inflammation Down-DGE",
                                  "TNF inflammation Down-DGE"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.58], polygon[GRID.polygon.59], polygon[GRID.polygon.60], polygon[GRID.polygon.61], text[GRID.text.62], text[GRID.text.63], lines[GRID.lines.64], text[GRID.text.65], lines[GRID.lines.66], text[GRID.text.67], text[GRID.text.68])

dev.off()

## quartz_off_screen 
##                 2

inf_overlap <- c(inf_overlap_up, inf_overlap_down)

It appears that 173/178 genes in total inflammation DGE overlap are dysregulated in the same direction between TCT and TNF, which is comforting. I will proceed with the remaining of the analysis using the total inflammation DGE list.

GO/KEGG/GSEA Analysis

Overlapping genes

GO analysis

target_gene <- inf_overlap
tct_all <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_I_V vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

tnf_all <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_I_V vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

detected_gene <- unique(c(tct_all$X1, tnf_all$X1))

# GO analysis
# Run GO enrichment analysis for biological process
ego_BP <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "BP", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_BP <- data.frame(ego_BP)

write.csv(cluster_summary_BP, "GO analysis/GO analysis_BP_overlap.csv")

# Run GO enrichment analysis for molecular function 
ego_MF <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "MF", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_MF <- data.frame(ego_MF)

write.csv(cluster_summary_MF, "GO analysis/GO analysis_MF_overlap.csv")

# Run GO enrichment analysis for cellular component 
ego_CC <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "CC", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_CC <- data.frame(ego_CC)

write.csv(cluster_summary_CC, "GO analysis/GO analysis_CC_overlap.csv")

Biological process

pdf("GO analysis/GO dotplot_BP_overlap.pdf",
           width = 8,
           height = 4)
dotplot(ego_BP, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_BP_overlap.pdf',
    width = 14,
    height = 14)
emapplot(ego_BP, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

Molecular Function

pdf('GO analysis/GO dotplot_MF_overlap.pdf',
    width = 9,
    height = 4)
dotplot(ego_MF, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_MF_overlap.pdf',
    width = 8,
    height = 8)
emapplot(ego_MF, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

Cellular component

pdf('GO analysis/GO dotplot_CC_overlap.pdf',
    width = 8,
    height = 4)
dotplot(ego_CC, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_CC_overlap.pdf',
    width = 8,
    height = 8)
emapplot(ego_CC, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

KEGG analysis

KEGG analysis needs Entrez ID.

# Convert Ensembl ID to Entrez ID
annotations_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = target_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

target_gene_entrez <- as.character(unique(annotations_entrez$ENTREZID[!is.na(annotations_entrez$ENTREZID)]))

annotations_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = detected_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

detected_gene_entrez <- as.character(unique(annotations_entrez$ENTREZID[!is.na(annotations_entrez$ENTREZID)]))

kk <- enrichKEGG(gene = target_gene_entrez, universe = detected_gene_entrez, organism = 'mmu')

## Reading KEGG annotation online:
## 
## Reading KEGG annotation online:

write.csv(as.data.frame(kk), "KEGG analysis/KEGG_overlap.csv")

dotplot(kk, showCategory = 50)

pdf("KEGG analysis/KEGG_overlap_inflammation_DGE.pdf",
    height = 6,
    width = 8)
dotplot(kk, showCategory = 20)
dev.off()

## quartz_off_screen 
##                 2

GSEA analysis

I only performed GSEA on the Hallmark geneset.

I am ranking the overlapping genes based on the average LFC between TCT and TNF model.

target_tct_LFC <- tct_inf$log2FoldChange[tct_inf$X1 %in% target_gene] 
target_tnf_LFC <- tnf_inf$log2FoldChange[tnf_inf$X1 %in% target_gene] 
target_LFC <- as.data.frame(cbind(target_gene, target_tct_LFC, target_tnf_LFC))

target_LFC$target_gene <- as.character(target_LFC$target_gene)
target_LFC$target_tct_LFC <- as.numeric(as.character(target_LFC$target_tct_LFC))
target_LFC$target_tnf_LFC <- as.numeric(as.character(target_LFC$target_tnf_LFC))

target_LFC$average_LFC <- rowMeans(target_LFC[,c(2:3)])

annotations_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = target_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

target_LFC <- inner_join(target_LFC, annotations_entrez, by = c("target_gene" = "GENEID"))
target_LFC <- target_LFC[!is.na(target_LFC$ENTREZID),]

load("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/Joint_Analysis/GSEA analysis/mouse_H_v5p2.rdata")
pathways <- Mm.H

gseadata <- target_LFC$ENTREZID
ranks <- target_LFC$average_LFC
names(ranks) <- target_LFC$ENTREZID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (0.65% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
write.csv(as.data.frame(fgseaRes), "GSEA analysis/GSEA_overlap.csv")
plotEnrichment(pathways[[fgseaRes$pathway[1]]], ranks) + labs(title = fgseaRes$pathway[1])

pdf('GSEA analysis/GSEA_overlapping_DGE.pdf',
    width = 8,
    height = 1)
plotGseaTable(pathways[fgseaRes$pathway], ranks, fgseaRes, gseaParam = 0.5) 
dev.off()

## quartz_off_screen 
##                 2

GSEA

TCT-specific genes

TCT specific DGEs

GO analysis

target_gene <- setdiff(tct_inf$X1, inf_overlap)
detected_gene <- tct_all$X1
# GO analysis
# Run GO enrichment analysis for biological process
ego_BP <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "BP", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_BP <- data.frame(ego_BP)

write.csv(cluster_summary_BP, "GO analysis/GO analysis_BP_TCT.csv")

# Run GO enrichment analysis for molecular function 
ego_MF <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "MF", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_MF <- data.frame(ego_MF)

write.csv(cluster_summary_MF, "GO analysis/GO analysis_MF_TCT.csv")

# Run GO enrichment analysis for cellular component 
ego_CC <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "CC", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_CC <- data.frame(ego_CC)

write.csv(cluster_summary_CC, "GO analysis/GO analysis_CC_TCT.csv")

Biological process

pdf('GO analysis/GO dotplot_BP_TCT.pdf',
    width = 10,
    height = 4)
dotplot(ego_BP, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_BP_TCT.pdf',
    width = 14,
    height = 14)
emapplot(ego_BP, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

Molecular Function

pdf('GO analysis/GO dotplot_MF_TCT.pdf',
    width = 8,
    height = 4)
dotplot(ego_MF, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_MF_TCT.pdf',
    width = 8,
    height = 8)
emapplot(ego_MF, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

Cellular component

pdf('GO analysis/GO dotplot_CC_TCT.pdf',
    width = 8,
    height = 4)
dotplot(ego_CC, showCategory = 10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_CC_TCT.pdf',
    width = 8,
    height = 8)
emapplot(ego_CC, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

KEGG analysis

KEGG analysis needs Entrez ID.

# Convert Ensembl ID to Entrez ID
annotations_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = target_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

target_gene_entrez <- as.character(unique(annotations_entrez$ENTREZID[!is.na(annotations_entrez$ENTREZID)]))

annotations_detected_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = detected_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

detected_gene_entrez <- as.character(unique(annotations_detected_entrez$ENTREZID[!is.na(annotations_detected_entrez$ENTREZID)]))

kk <- enrichKEGG(gene = target_gene_entrez, universe = detected_gene_entrez, organism = 'mmu')
write.csv(as.data.frame(kk), "KEGG analysis/KEGG_TCT.csv")
dotplot(kk, showCategory = 20)

pdf("KEGG analysis/TCT_specific_KEGG analysis.pdf",
    height = 4,
    width = 10)
dotplot(kk, showCategory = 20)
dev.off()

## quartz_off_screen 
##                 2

GSEA analysis

I only performed GSEA on the Hallmark geneset.

target_LFC <- tct_inf[tct_inf$X1 %in% target_gene,]

annotations_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = target_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

target_LFC <- inner_join(target_LFC, annotations_entrez, by = c("X1" = "GENEID"))
target_LFC <- target_LFC[!is.na(target_LFC$ENTREZID),]

gseadata <- target_LFC$ENTREZID
ranks <- target_LFC$log2FoldChange
names(ranks) <- target_LFC$ENTREZID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (2.52% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

# number of GSEA pathway enriched
dim(fgseaRes)[1] 

## [1] 0

TCT exluding overlapping inflammation DGEs

GSEA analysis

I only performed GSEA on the Hallmark geneset.

The ranking is done usin the default signal to noise ratio from GSEA software.

TCT_iv_niv_norm <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/normalized_raw_gene_counts_I_V vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_6 = col_double(),
##   NI_V_7 = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_V_6 = col_double()
## )

target_norm <- TCT_iv_niv_norm[!TCT_iv_niv_norm$X1 %in% inf_overlap, ]
write_csv(target_norm, "GSEA analysis/TCT excluding overlap/normalized_counts.csv")

tct_rank <- read_csv("GSEA analysis/TCT excluding overlap/Hallmark/TCT_rank.csv")

## Parsed with column specification:
## cols(
##   NAME = col_character(),
##   DESCRIPTION = col_character(),
##   GENE_SYMBOL = col_character(),
##   GENE_TITLE = col_character(),
##   SCORE = col_double()
## )

human2mouse <- read_csv("Human_2_mouse_symbol.csv")

## Parsed with column specification:
## cols(
##   `Human Marker Symbol` = col_character(),
##   `Human Entrez Gene ID` = col_double(),
##   `Mouse Marker Symbol` = col_character(),
##   `MGI Marker Accession ID` = col_character()
## )

mouseID <- read_csv("mouse_symbol_Ensembl.csv")

## Parsed with column specification:
## cols(
##   `MGI Accession ID` = col_character(),
##   `Marker Symbol` = col_character(),
##   `EntrezGene ID` = col_character(),
##   `Ensembl Gene ID` = col_character(),
##   `HomoloGene ID` = col_character()
## )

tct_rank <- inner_join(tct_rank, human2mouse[,c(1,3)], by = c("NAME" = "Human Marker Symbol"))
tct_rank <- inner_join(tct_rank, mouseID[,c(2,3)], by = c("Mouse Marker Symbol" = "Marker Symbol"))

target_LFC <- tct_rank[!is.na(tct_rank$`EntrezGene ID`),]

gseadata <- target_LFC$`EntrezGene ID`
ranks <- target_LFC$SCORE
names(ranks) <- target_LFC$`EntrezGene ID`
ranks_tct <- sort(ranks, decreasing = T)
fgseaRes_tct <- fgsea(pathways, ranks_tct, minSize=15, maxSize = 500)

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (1.95% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.

## Warning in fgseaMultilevel(...): For some pathways, in reality P-values are less
## than 1e-10. You can set the `eps` argument to zero for better estimation.

fgseaRes_tct_sig <- fgseaRes_tct[fgseaRes_tct$pval < 0.05 & fgseaRes_tct$padj < 0.1, ]

fgseaRes_tct_sig <- fgseaRes_tct_sig[order(fgseaRes_tct_sig$NES, decreasing = TRUE), ]

# number of GSEA pathway enriched
dim(fgseaRes_tct_sig)[1]

## [1] 28

suppressMessages(fgseaRes_tct[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes_tct$leadingEdge <- as.character(fgseaRes_tct$leadingEdge)
write.csv(as.data.frame(fgseaRes_tct), "GSEA analysis/GSEA_TCT_excluding overlap.csv")

pdf('GSEA analysis/GSEA_TCT_excluding overlap.pdf',
    width = 12,
    height = 10)
plotGseaTable(pathways[fgseaRes_tct_sig$pathway], ranks_tct, fgseaRes_tct_sig, gseaParam = 0.5)
dev.off()

## quartz_off_screen 
##                 2

GSEA

TNF-specific DGE

GO analysis

target_gene <- setdiff(tnf_inf$X1, inf_overlap)
detected_gene <- tnf_all$X1

# GO analysis
# Run GO enrichment analysis for biological process
ego_BP <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "BP", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_BP <- data.frame(ego_BP)

write.csv(cluster_summary_BP, "GO analysis/GO analysis_BP_TNF.csv")

# Run GO enrichment analysis for molecular function 
ego_MF <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "MF", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_MF <- data.frame(ego_MF)

write.csv(cluster_summary_MF, "GO analysis/GO analysis_MF_TNF.csv")

# Run GO enrichment analysis for cellular component 
ego_CC <- enrichGO(gene = target_gene, 
                universe = detected_gene,
                keyType = "ENSEMBL",
                OrgDb = org.Mm.eg.db, 
                ont = "CC", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                readable = TRUE)

# Output results from GO analysis to a table
cluster_summary_CC <- data.frame(ego_CC)

write.csv(cluster_summary_CC, "GO analysis/GO analysis_CC_TNF.csv")

Biological process

pdf('GO analysis/GO dotplot_BP_TNF.pdf',
    width = 8,
    height = 4)
dotplot(ego_BP, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_BP_TNF.pdf',
    width = 14,
    height = 14)
emapplot(ego_BP, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

Molecular Function

pdf('GO analysis/GO dotplot_MF_TNF.pdf',
    width = 8,
    height = 4)
dotplot(ego_MF, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_MF_TNF.pdf',
    width = 14,
    height = 14)
emapplot(ego_MF, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

Cellular component

pdf('GO analysis/GO dotplot_CC_TNF.pdf',
    width = 8,
    height = 4)
dotplot(ego_CC, showCategory=10)
dev.off()

## quartz_off_screen 
##                 2

pdf('GO analysis/GO enrichment_CC_TNF.pdf',
    width = 14,
    height = 14)
emapplot(ego_CC, showCategory = 50)
dev.off()

## quartz_off_screen 
##                 2

Final output is following:

KEGG analysis

KEGG analysis needs Entrez ID.

# Convert Ensembl ID to Entrez ID
annotations_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = target_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

target_gene_entrez <- as.character(unique(annotations_entrez$ENTREZID[!is.na(annotations_entrez$ENTREZID)]))

annotations_detected_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = detected_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

detected_gene_entrez <- as.character(unique(annotations_detected_entrez$ENTREZID[!is.na(annotations_detected_entrez$ENTREZID)]))

kk <- enrichKEGG(gene = target_gene_entrez, universe = detected_gene_entrez, organism = 'mmu')
write.csv(as.data.frame(kk), "KEGG analysis/KEGG_TCT.csv")
dotplot(kk, showCategory = 20)

pdf("KEGG analysis/TNF_specific_KEGG analysis.pdf",
    height = 6,
    width = 10)
dotplot(kk, showCategory = 20)
dev.off()

## quartz_off_screen 
##                 2

GSEA analysis

I only performed GSEA on the Hallmark geneset.

target_LFC <- tnf_inf[tnf_inf$X1 %in% target_gene,]

annotations_entrez <- AnnotationDbi::select(EnsDb.Mmusculus.v79,
                                           keys = target_gene,
                                           columns = c("ENTREZID"),
                                           keytype = "GENEID")

target_LFC <- inner_join(target_LFC, annotations_entrez, by = c("X1" = "GENEID"))
target_LFC <- target_LFC[!is.na(target_LFC$ENTREZID),]

gseadata <- target_LFC$ENTREZID
ranks <- target_LFC$log2FoldChange
names(ranks) <- target_LFC$ENTREZID
ranks <- sort(ranks, decreasing = T)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500)

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (3.27% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.

fgseaRes <- fgseaRes[fgseaRes$pval< 0.05 & fgseaRes$padj < 0.1,]
fgseaRes <- fgseaRes[order(fgseaRes$NES, decreasing = TRUE), ]

# number of GSEA pathway enriched
dim(fgseaRes)[1]

## [1] 13

suppressMessages(fgseaRes[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes$leadingEdge <- as.character(fgseaRes$leadingEdge)
write.csv(as.data.frame(fgseaRes), "GSEA analysis/GSEA_TNF_specific.csv")
pdf('GSEA analysis/GSEA_TNF_specific_DGE.pdf',
    width = 10,
    height = 5)
plotGseaTable(pathways[fgseaRes$pathway], ranks, fgseaRes, gseaParam = 0.5)
dev.off()

## quartz_off_screen 
##                 2

GSEA

TNF exluding overlapping inflammation DGEs

GSEA analysis

I only performed GSEA on the Hallmark geneset.

The ranking is done using the default signal to noise ratio from GSEA software.

TNF_iv_niv_norm <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/normalized_raw_gene_counts_I_V vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_3 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_5 = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double()
## )

target_norm <- TNF_iv_niv_norm[!TNF_iv_niv_norm$X1 %in% inf_overlap, ]
write_csv(target_norm, "GSEA analysis/TNF excluding overlap/normalized_counts.csv")

tnf_rank <- read_csv("GSEA analysis/TNF excluding overlap/Hallmark/TNF_rank.csv")

## Parsed with column specification:
## cols(
##   NAME = col_character(),
##   DESCRIPTION = col_character(),
##   GENE_SYMBOL = col_character(),
##   GENE_TITLE = col_character(),
##   SCORE = col_double()
## )

human2mouse <- read_csv("Human_2_mouse_symbol.csv")

## Parsed with column specification:
## cols(
##   `Human Marker Symbol` = col_character(),
##   `Human Entrez Gene ID` = col_double(),
##   `Mouse Marker Symbol` = col_character(),
##   `MGI Marker Accession ID` = col_character()
## )

mouseID <- read_csv("mouse_symbol_Ensembl.csv")

## Parsed with column specification:
## cols(
##   `MGI Accession ID` = col_character(),
##   `Marker Symbol` = col_character(),
##   `EntrezGene ID` = col_character(),
##   `Ensembl Gene ID` = col_character(),
##   `HomoloGene ID` = col_character()
## )

tnf_rank <- inner_join(tnf_rank, human2mouse[,c(1,3)], by = c("NAME" = "Human Marker Symbol"))
tnf_rank <- inner_join(tnf_rank, mouseID[,c(2,3)], by = c("Mouse Marker Symbol" = "Marker Symbol"))

target_LFC <- tnf_rank[!is.na(tnf_rank$`EntrezGene ID`),]

gseadata <- target_LFC$`EntrezGene ID`
ranks <- target_LFC$SCORE
names(ranks) <- target_LFC$`EntrezGene ID`
ranks_tnf <- sort(ranks, decreasing = T)
fgseaRes_tnf <- fgsea(pathways, ranks_tnf, minSize=15, maxSize = 500)

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (1.99% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There are duplicate gene names, fgsea may produce unexpected
## results.

## Warning in fgseaMultilevel(...): For some pathways, in reality P-values are less
## than 1e-10. You can set the `eps` argument to zero for better estimation.

fgseaRes_tnf_sig <- fgseaRes_tnf[fgseaRes_tnf$pval < 0.05 & fgseaRes_tnf$padj < 0.1, ]

fgseaRes_tnf_sig <- fgseaRes_tnf_sig[order(fgseaRes_tnf_sig$NES, decreasing = TRUE), ]

# number of GSEA pathway enriched
dim(fgseaRes_tnf_sig)[1]

## [1] 25

suppressMessages(fgseaRes_tnf[, leadingEdge := lapply(leadingEdge, mapIds, x=org.Mm.eg.db, keytype="ENTREZID", column="SYMBOL")])
fgseaRes_tnf$leadingEdge <- as.character(fgseaRes_tnf$leadingEdge)
write.csv(as.data.frame(fgseaRes_tnf), "GSEA analysis/GSEA_TNF_excluding overlap.csv")

pdf('GSEA analysis/GSEA_TNF_excluding overlap.pdf',
    width = 12,
    height = 10)
plotGseaTable(pathways[fgseaRes_tnf_sig$pathway], ranks_tct, fgseaRes_tnf_sig, gseaParam = 0.5)
dev.off()

## quartz_off_screen 
##                 2

GSEA

Joint analysis of TCT-specific and TNF-specific GSEA

This analysis is done using RNA-Seq excluding overlapping inflammation related DGE.

I would like to look at what is the enriched GSEA hallmark pathways that are unqiue to TCT or TNF inflammation and the overlap.

TCT_specific_gsea <- setdiff(fgseaRes_tct_sig$pathway, fgseaRes_tnf_sig$pathway)
TNF_specific_gsea <- setdiff(fgseaRes_tnf_sig$pathway, fgseaRes_tct_sig$pathway)
overlap_gsea <- intersect(fgseaRes_tct_sig$pathway, fgseaRes_tnf_sig$pathway)


# Now I need to check if all the overlapping GSEA has the same direction of enrichment
for (i in overlap_gsea) {
  tct_nes <- fgseaRes_tct_sig$NES[fgseaRes_tct_sig$pathway == i]
  tnf_nes <- fgseaRes_tnf_sig$NES[fgseaRes_tnf_sig$pathway == i]
  dir <- tct_nes * tnf_nes
  if (dir < 0) {
    print(i)
  }
}

## [1] "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"
## [1] "HALLMARK_APICAL_JUNCTION"

# plot out the overlapping pathways
pdf('GSEA analysis/Overlapping_GSEA_TCT.pdf',
    width = 12,
    height = 6)
plotGseaTable(pathways[overlap_gsea], ranks_tct, fgseaRes_tct_sig, gseaParam = 0.5) +
dev.off()

## integer(0)

GSEA

pdf('GSEA analysis/Overlapping_GSEA_TNF.pdf',
    width = 12,
    height = 6)
plotGseaTable(pathways[overlap_gsea], ranks_tnf, fgseaRes_tnf_sig, gseaParam = 0.5)
dev.off()

## quartz_off_screen 
##                 2

GSEA

# plot out the model-specific pathways
## TCT specific
pdf('GSEA analysis/GSEA_TCT_specific.pdf',
    width = 12,
    height = 4)
plotGseaTable(pathways[TCT_specific_gsea], ranks_tct, fgseaRes_tct_sig, gseaParam = 0.5)
dev.off()

## quartz_off_screen 
##                 2

GSEA

# plot out the model-specific pathways
## TNF specific
pdf('GSEA analysis/GSEA_TNF_specific.pdf',
    width = 12,
    height = 4)
plotGseaTable(pathways[TNF_specific_gsea], ranks_tnf, fgseaRes_tnf_sig, gseaParam = 0.5)
dev.off()

## quartz_off_screen 
##                 2

GSEA

Identify MK2 target genes

There are 4 treatments of MK2i that can be used for comparison. TCT and TNF model each have MK2i treatment in inflammaed and non-inflammed conditions.

Again I am using 0.1 as the adjusted p-value cutoff.

First let's examine the TCT model.

# Non-inflammed condition
tct_ni_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Significant genes_NI_MKI vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_6 = col_double(),
##   NI_V_7 = col_double(),
##   NI_MKI_1 = col_double(),
##   NI_MKI_2 = col_double(),
##   NI_MKI_3 = col_double(),
##   NI_MKI_4 = col_double(),
##   NI_MKI_5 = col_double(),
##   NI_MKI_6 = col_double()
## )

tct_ni_MK2i_DGE <- tct_ni_MK2i$X1

# Inflammed condition
tct_i_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Significant genes_I_MKI vs I_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_V_6 = col_double(),
##   I_MKI_1 = col_double(),
##   I_MKI_2 = col_double(),
##   I_MKI_3 = col_double(),
##   I_MKI_4 = col_double(),
##   I_MKI_5 = col_double(),
##   I_MKI_6 = col_double()
## )

tct_i_MK2i_DGE <- tct_i_MK2i$X1

# combined condition
tct_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Significant genes_MKI vs V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   X1 = col_character()
## )

## See spec(...) for full column specifications.

tct_MK2i_DGE <- tct_MK2i$X1

# Examine the overlap between NI and I conditions and also combined analysis

n12 <- intersect(tct_ni_MK2i_DGE, tct_i_MK2i_DGE)
n23 <- intersect(tct_i_MK2i_DGE, tct_MK2i_DGE)
n13 <- intersect(tct_ni_MK2i_DGE, tct_MK2i_DGE)
n123 <- intersect(n12, tct_MK2i_DGE)

grid.newpage()
draw.triple.venn(area1 = length(tct_ni_MK2i_DGE),
                 area2 = length(tct_i_MK2i_DGE),
                 area3 = length(tct_MK2i_DGE),
                 n12 = length(n12),
                 n23 = length(n23),
                 n13 = length(n13),
                 n123 = length(n123),
                 catergory <- c("TCT MK2i target\nnon-inflammation",
                                "TCT MK2i target\ninflammation",
                                "TCT MK2i target\ncombined"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5568], polygon[GRID.polygon.5569], polygon[GRID.polygon.5570], polygon[GRID.polygon.5571], polygon[GRID.polygon.5572], polygon[GRID.polygon.5573], text[GRID.text.5574], text[GRID.text.5575], text[GRID.text.5576], text[GRID.text.5577], text[GRID.text.5578], text[GRID.text.5579], text[GRID.text.5580])

pdf("TCT-i_vs_ni_MK2i target.pdf",
    height = 8,
    width = 12)
grid.newpage()
draw.triple.venn(area1 = length(tct_ni_MK2i_DGE),
                 area2 = length(tct_i_MK2i_DGE),
                 area3 = length(tct_MK2i_DGE),
                 n12 = length(n12),
                 n23 = length(n23),
                 n13 = length(n13),
                 n123 = length(n123),
                 catergory <- c("TCT MK2i target\nnon-inflammation",
                                "TCT MK2i target\ninflammation",
                                "TCT MK2i target\ncombined"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5581], polygon[GRID.polygon.5582], polygon[GRID.polygon.5583], polygon[GRID.polygon.5584], polygon[GRID.polygon.5585], polygon[GRID.polygon.5586], text[GRID.text.5587], text[GRID.text.5588], text[GRID.text.5589], text[GRID.text.5590], text[GRID.text.5591], text[GRID.text.5592], text[GRID.text.5593])

dev.off()

## quartz_off_screen 
##                 2

It appears that there is no overlapping MK2i targets between the inflammed and non-inflammed TCT models.

Now let's examine the TNF model.

# Non-inflammed condition
tnf_ni_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Significant genes_NI_MKI vs NI_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_3 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_5 = col_double(),
##   NI_MKI_1 = col_double(),
##   NI_MKI_2 = col_double(),
##   NI_MKI_3 = col_double(),
##   NI_MKI_4 = col_double(),
##   NI_MKI_5 = col_double()
## )

tnf_ni_MK2i_DGE <- tnf_ni_MK2i$X1

# Inflammed condition
tnf_i_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Significant genes_I_MKI vs I_V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_MKI_1 = col_double(),
##   I_MKI_2 = col_double(),
##   I_MKI_3 = col_double(),
##   I_MKI_4 = col_double(),
##   I_MKI_5 = col_double(),
##   I_MKI_6 = col_double()
## )

tnf_i_MK2i_DGE <- tnf_i_MK2i$X1

# combined condition
tnf_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Significant genes_MKI vs V.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   X1 = col_character()
## )

## See spec(...) for full column specifications.

tnf_MK2i_DGE <- tnf_MK2i$X1

# Examine the overlap between NI and I conditions
n12 <- intersect(tnf_ni_MK2i_DGE, tnf_i_MK2i_DGE)
n23 <- intersect(tnf_i_MK2i_DGE, tnf_MK2i_DGE)
n13 <- intersect(tnf_ni_MK2i_DGE, tnf_MK2i_DGE)
n123 <- intersect(n12, tnf_MK2i_DGE)

grid.newpage()
draw.triple.venn(area1 = length(tnf_ni_MK2i_DGE),
                 area2 = length(tnf_i_MK2i_DGE),
                 area3 = length(tnf_MK2i_DGE),
                 n12 = length(n12),
                 n23 = length(n23),
                 n13 = length(n13),
                 n123 = length(n123),
                 catergory <- c("TNF MK2i target\nnon-inflammation",
                                "TNF MK2i target\ninflammation",
                                "TNF MK2i target\ncombined"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5594], polygon[GRID.polygon.5595], polygon[GRID.polygon.5596], polygon[GRID.polygon.5597], polygon[GRID.polygon.5598], polygon[GRID.polygon.5599], text[GRID.text.5600], text[GRID.text.5601], text[GRID.text.5602], text[GRID.text.5603], text[GRID.text.5604], text[GRID.text.5605], text[GRID.text.5606], text[GRID.text.5607])

pdf("TNF-i_vs_ni_MK2i target.pdf",
    height = 8,
    width = 12)
grid.newpage()
draw.triple.venn(area1 = length(tnf_ni_MK2i_DGE),
                 area2 = length(tnf_i_MK2i_DGE),
                 area3 = length(tnf_MK2i_DGE),
                 n12 = length(n12),
                 n23 = length(n23),
                 n13 = length(n13),
                 n123 = length(n123),
                 catergory <- c("TNF MK2i target\nnon-inflammation",
                                "TNF MK2i target\ninflammation",
                                "TNF MK2i target\ncombined"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5608], polygon[GRID.polygon.5609], polygon[GRID.polygon.5610], polygon[GRID.polygon.5611], polygon[GRID.polygon.5612], polygon[GRID.polygon.5613], text[GRID.text.5614], text[GRID.text.5615], text[GRID.text.5616], text[GRID.text.5617], text[GRID.text.5618], text[GRID.text.5619], text[GRID.text.5620], text[GRID.text.5621])

dev.off()

## quartz_off_screen 
##                 2

Now let's take a looks and see if there is any overlap between TCT and TNF models.

# Inflammed conditions
i_MK2i_overlap <- intersect(tnf_i_MK2i_DGE, tct_i_MK2i_DGE)
grid.newpage()
draw.pairwise.venn(length(tnf_i_MK2i_DGE),
                   length(tct_i_MK2i_DGE),
                   length(i_MK2i_overlap),
                   catergory <- c("MK2i target in\ninflammed TNF",
                                  "MK2i target in\ninflammed TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(135, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5622], polygon[GRID.polygon.5623], polygon[GRID.polygon.5624], polygon[GRID.polygon.5625], text[GRID.text.5626], text[GRID.text.5627], text[GRID.text.5628], text[GRID.text.5629])

pdf("Inflammed_TCT_vs_TNF_MKi_target.pdf",
    height = 8,
    width = 10)
draw.pairwise.venn(length(tnf_i_MK2i_DGE),
                   length(tct_i_MK2i_DGE),
                   length(i_MK2i_overlap),
                   catergory <- c("MK2i target in\ninflammed TNF",
                                  "MK2i target in\ninflammed TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(135, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5630], polygon[GRID.polygon.5631], polygon[GRID.polygon.5632], polygon[GRID.polygon.5633], text[GRID.text.5634], text[GRID.text.5635], text[GRID.text.5636], text[GRID.text.5637])

dev.off()

## quartz_off_screen 
##                 2

# Non-inflammed conditions
ni_MK2i_overlap <- intersect(tnf_ni_MK2i_DGE, tct_ni_MK2i_DGE)
grid.newpage()
draw.pairwise.venn(length(tnf_ni_MK2i_DGE),
                   length(tct_ni_MK2i_DGE),
                   length(ni_MK2i_overlap),
                   catergory <- c("MK2i target in\nnon-inflammed TNF",
                                  "MK2i target in\nnon-inflammed TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5638], polygon[GRID.polygon.5639], polygon[GRID.polygon.5640], polygon[GRID.polygon.5641], text[GRID.text.5642], text[GRID.text.5643], text[GRID.text.5644], text[GRID.text.5645])

pdf("Non-inflammed_TCT_vs_TNF_MKi_target.pdf",
    height = 8,
    width = 10)
draw.pairwise.venn(length(tnf_ni_MK2i_DGE),
                   length(tct_ni_MK2i_DGE),
                   length(ni_MK2i_overlap),
                   catergory <- c("MK2i target in\nnon-inflammed TNF",
                                  "MK2i target in\nnon-inflammed TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5646], polygon[GRID.polygon.5647], polygon[GRID.polygon.5648], polygon[GRID.polygon.5649], text[GRID.text.5650], text[GRID.text.5651], text[GRID.text.5652], text[GRID.text.5653])

dev.off()

## quartz_off_screen 
##                 2

# Combined condition
MK2i_overlap <- intersect(tnf_MK2i_DGE, tct_MK2i_DGE)
grid.newpage()
draw.pairwise.venn(length(tnf_MK2i_DGE),
                   length(tct_MK2i_DGE),
                   length(MK2i_overlap),
                   catergory <- c("MK2i target in\ncombined TNF",
                                  "MK2i target in\ncombined TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(135, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5654], polygon[GRID.polygon.5655], polygon[GRID.polygon.5656], polygon[GRID.polygon.5657], text[GRID.text.5658], text[GRID.text.5659], text[GRID.text.5660], text[GRID.text.5661])

So it appears that none of the 4 sets of DGEs have any overlap with each other. That is a little disturbing.

Identify Biomarkers of effective response to MK2i.

Biomarkers are defined as genes that are targets of MK2i and are dysregulated in only TCT model of inflammation but not TNF model of inflammation.

First we will be more stringent, I will only use targets of MK2i in inflammed TCT model.

tct_biomarker_raw <- intersect(tct_i_MK2i_DGE, tct_inf_DGE)

tct_biomarker <- setdiff(tct_biomarker_raw, tnf_inf_DGE)

n23 <- intersect(tct_i_MK2i_DGE, tnf_inf_DGE)

n123 <- intersect(inf_overlap, tct_i_MK2i_DGE)

grid.newpage()
draw.triple.venn(area1 = length(tct_inf_DGE),
                 area2 = length(tnf_inf_DGE),
                 area3 = length(tct_i_MK2i_DGE),
                 n12 = length(inf_overlap),
                 n23 = length(n23),
                 n13 = length(tct_biomarker_raw),
                 n123 = length(n123),
                 catergory <- c("TCT inflammation DGE",
                                "TNF inflammation DGE",
                                "MK2i target in inflammed TCT"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5662], polygon[GRID.polygon.5663], polygon[GRID.polygon.5664], polygon[GRID.polygon.5665], polygon[GRID.polygon.5666], polygon[GRID.polygon.5667], text[GRID.text.5668], text[GRID.text.5669], text[GRID.text.5670], text[GRID.text.5671], text[GRID.text.5672], text[GRID.text.5673], text[GRID.text.5674], text[GRID.text.5675])

pdf("MK2i responsiveness marker.pdf",
    height = 8,
    width = 10)
draw.triple.venn(area1 = length(tct_inf_DGE),
                 area2 = length(tnf_inf_DGE),
                 area3 = length(tct_i_MK2i_DGE),
                 n12 = length(inf_overlap),
                 n23 = length(n23),
                 n13 = length(tct_biomarker_raw),
                 n123 = length(n123),
                 catergory <- c("TCT inflammation DGE",
                                "TNF inflammation DGE",
                                "MK2i target in inflammed TCT"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5676], polygon[GRID.polygon.5677], polygon[GRID.polygon.5678], polygon[GRID.polygon.5679], polygon[GRID.polygon.5680], polygon[GRID.polygon.5681], text[GRID.text.5682], text[GRID.text.5683], text[GRID.text.5684], text[GRID.text.5685], text[GRID.text.5686], text[GRID.text.5687], text[GRID.text.5688], text[GRID.text.5689])

dev.off()

## quartz_off_screen 
##                 2

So there are three genes that popped out as potential biomarkers. List them here:

tct_biomarker

## character(0)

We want to make sure the LFC are the opposite for these three genes in TCT inflammation DGE list and MK2i target in inflammed TCT.

tct_inf_DGE_LFC <- tct_inf$log2FoldChange[tct_inf$X1 %in% tct_biomarker]
tct_i_MK2i_LFC <- tct_i_MK2i$log2FoldChange[tct_i_MK2i$X1 %in% tct_biomarker]
biomarker_table <- cbind(tct_biomarker, tct_inf_DGE_LFC, tct_i_MK2i_LFC)
colnames(biomarker_table) <- c("Ensembl_ID", "LFC(TCT_inflammation_related)", "LFC(TCT_inflammed_MK2i_target)")
biomarker_table

##      Ensembl_ID LFC(TCT_inflammation_related) LFC(TCT_inflammed_MK2i_target)

Now I will combined all 4 sets of MK2i targets and use that as MK2i targets.

mki_target <- c(tnf_i_MK2i_DGE, tnf_ni_MK2i_DGE, tct_i_MK2i_DGE, tct_ni_MK2i_DGE)

tct_biomarker_raw <- intersect(mki_target, tct_inf_DGE)

tct_biomarker <- setdiff(tct_biomarker_raw, tnf_inf_DGE)

n23 <- intersect(mki_target, tnf_inf_DGE)

n123 <- intersect(inf_overlap, mki_target)

grid.newpage()
draw.triple.venn(area1 = length(tct_inf_DGE),
                 area2 = length(tnf_inf_DGE),
                 area3 = length(mki_target),
                 n12 = length(inf_overlap),
                 n23 = length(n23),
                 n13 = length(tct_biomarker_raw),
                 n123 = length(n123),
                 catergory <- c("TCT inflammation DGE",
                                "TNF inflammation DGE",
                                "All MK2i targets"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5690], polygon[GRID.polygon.5691], polygon[GRID.polygon.5692], polygon[GRID.polygon.5693], polygon[GRID.polygon.5694], polygon[GRID.polygon.5695], text[GRID.text.5696], text[GRID.text.5697], text[GRID.text.5698], text[GRID.text.5699], text[GRID.text.5700], text[GRID.text.5701], text[GRID.text.5702], text[GRID.text.5703], text[GRID.text.5704])

So there are three genes that popped out as potential biomarkers. List them here:

tct_biomarker

## character(0)

We want to make sure the LFC are the opposite for these three genes in TCT inflammation DGE list and MK2i target in inflammed TCT.

tct_inf_DGE_LFC <- tct_inf$log2FoldChange[tct_inf$X1 %in% tct_biomarker]
tct_i_MK2i_LFC <- tct_i_MK2i$log2FoldChange[tct_i_MK2i$X1 %in% tct_biomarker]
biomarker_table <- cbind(tct_biomarker, tct_inf_DGE_LFC, tct_i_MK2i_LFC)
colnames(biomarker_table) <- c("Ensembl_ID", "LFC(TCT_inflammation_related)", "LFC(TCT_inflammed_MK2i_target)")
biomarker_table

##      Ensembl_ID LFC(TCT_inflammation_related) LFC(TCT_inflammed_MK2i_target)

It appears that MK2i are downregulating these three genes, while inflammation upregulates them in the TCT model. This indicates that they could be the genes that mediate the responsiveness to MK2i treatment.

Identify Biomarkers of effective response to MK2i using filtered data (analysis without I_V3 sample)

tct_inf_DGE_filtered_data <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Significant genes_I_V vs NI_V_filtered.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_6 = col_double(),
##   NI_V_7 = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_V_6 = col_double()
## )

tct_inf_DGE_filtered <- tct_inf_DGE_filtered_data$X1

tct_i_MK2i_DGE_filtered_data <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Significant genes_I_MKI vs I_V_filtered.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_V_6 = col_double(),
##   I_MKI_1 = col_double(),
##   I_MKI_2 = col_double(),
##   I_MKI_3 = col_double(),
##   I_MKI_4 = col_double(),
##   I_MKI_5 = col_double(),
##   I_MKI_6 = col_double()
## )

tct_i_MK2i_DGE_filtered <- tct_i_MK2i_DGE_filtered_data$X1

tct_biomarker_raw <- intersect(tct_i_MK2i_DGE_filtered, tct_inf_DGE_filtered)

tct_biomarker <- setdiff(tct_biomarker_raw, tnf_inf_DGE)

inf_overlap_filtered <- intersect(tct_inf_DGE_filtered,tnf_inf_DGE)

n23 <- intersect(tct_i_MK2i_DGE_filtered, tnf_inf_DGE)

n123 <- intersect(inf_overlap, tct_i_MK2i_DGE_filtered)

grid.newpage()
draw.triple.venn(area1 = length(tct_inf_DGE_filtered),
                 area2 = length(tnf_inf_DGE),
                 area3 = length(tct_i_MK2i_DGE_filtered),
                 n12 = length(inf_overlap),
                 n23 = length(n23),
                 n13 = length(tct_biomarker_raw),
                 n123 = length(n123),
                 catergory <- c("TCT inflammation DGE",
                                "TNF inflammation DGE",
                                "MK2i target in inflammed TCT"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5705], polygon[GRID.polygon.5706], polygon[GRID.polygon.5707], polygon[GRID.polygon.5708], polygon[GRID.polygon.5709], polygon[GRID.polygon.5710], text[GRID.text.5711], text[GRID.text.5712], text[GRID.text.5713], text[GRID.text.5714], text[GRID.text.5715], text[GRID.text.5716], text[GRID.text.5717], text[GRID.text.5718], text[GRID.text.5719])

So there are three genes that popped out as potential biomarkers. List them here:

tct_biomarker

## [1] "ENSMUSG00000039199"

We want to make sure the LFC are the opposite for these three genes in TCT inflammation DGE list and MK2i target in inflammed TCT.

tct_inf_DGE_LFC <- tct_inf_DGE_filtered_data$log2FoldChange[tct_inf_DGE_filtered_data$X1 %in% tct_biomarker]
tct_i_MK2i_LFC <- tct_i_MK2i_DGE_filtered_data$log2FoldChange[tct_i_MK2i_DGE_filtered_data$X1 %in% tct_biomarker]
biomarker_table <- cbind(tct_biomarker, tct_inf_DGE_LFC, tct_i_MK2i_LFC)
colnames(biomarker_table) <- c("Ensembl_ID", "LFC(TCT_inflammation_related)", "LFC(TCT_inflammed_MK2i_target)")
biomarker_table

##      Ensembl_ID           LFC(TCT_inflammation_related)
## [1,] "ENSMUSG00000039199" "-4.77295995852485"          
##      LFC(TCT_inflammed_MK2i_target)
## [1,] "0.720290930939497"

Transcript level analysis

I am repeating the same analysis using the transcript level analysis result. ## Identify inflammation related DGE First I need to identify DEG in inflammation in TCT and TNF models.

# DGE in inflammation using I_V and NI_V comparison
tct_inf <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_I_V vs NI_V_transcript_level.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

tct_inf_DGE <- tct_inf$X1[tct_inf$padj < 0.05]
tct_inf_DGE_up <- tct_inf$X1[tct_inf$padj < 0.05 & tct_inf$log2FoldChange > 0]
tct_inf_DGE_down <- tct_inf$X1[tct_inf$padj < 0.05 & tct_inf$log2FoldChange < 0]

tnf_inf <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_I_V vs NI_V_transcript_level.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

tnf_inf_DGE <- tnf_inf$X1[tnf_inf$padj < 0.05]
tnf_inf_DGE_up <- tnf_inf$X1[tnf_inf$padj < 0.05 & tnf_inf$log2FoldChange > 0]
tnf_inf_DGE_down <- tnf_inf$X1[tnf_inf$padj < 0.05 & tnf_inf$log2FoldChange < 0]

# Find out how big is the overlap
## Total overlap
inf_overlap <-intersect(tct_inf_DGE, tnf_inf_DGE)

grid.newpage()
draw.pairwise.venn(length(tct_inf_DGE),
                   length(tnf_inf_DGE),
                   length(inf_overlap),
                   catergory <- c("TCT inflammation DGE_transcript",
                                  "TNF inflammation DGE_transcript"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5720], polygon[GRID.polygon.5721], polygon[GRID.polygon.5722], polygon[GRID.polygon.5723], text[GRID.text.5724], text[GRID.text.5725], text[GRID.text.5726], lines[GRID.lines.5727], text[GRID.text.5728], text[GRID.text.5729])

## Overlap in up-regulated genes
inf_overlap_up <-intersect(tct_inf_DGE_up, tnf_inf_DGE_up)

grid.newpage()
draw.pairwise.venn(length(tct_inf_DGE_up),
                   length(tnf_inf_DGE_up),
                   length(inf_overlap_up),
                   catergory <- c("TCT inflammation Up-DGE_transcript",
                                  "TNF inflammation Up-DGE_transcript"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5730], polygon[GRID.polygon.5731], polygon[GRID.polygon.5732], polygon[GRID.polygon.5733], text[GRID.text.5734], text[GRID.text.5735], text[GRID.text.5736], lines[GRID.lines.5737], text[GRID.text.5738], text[GRID.text.5739])

## Overlap in down-regulated genes
inf_overlap_down <-intersect(tct_inf_DGE_down, tnf_inf_DGE_down)

grid.newpage()
draw.pairwise.venn(length(tct_inf_DGE_down),
                   length(tnf_inf_DGE_down),
                   length(inf_overlap_down),
                   catergory <- c("TCT inflammation Down-DGE_transcript",
                                  "TNF inflammation Down-DGE_transcript"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5740], polygon[GRID.polygon.5741], polygon[GRID.polygon.5742], polygon[GRID.polygon.5743], text[GRID.text.5744], text[GRID.text.5745], text[GRID.text.5746], lines[GRID.lines.5747], text[GRID.text.5748], text[GRID.text.5749])

It appears that 82/90 transcripts in total inflammation DGE overlap are dysregulated in the same direction between TCT and TNF, which is comforting. I will proceed with the remaining of the analysis using the total inflammation DGE list.

Identify MK2 target genes

There are 4 treatments of MK2i that can be used for comparison. TCT and TNF model each have MK2i treatment in inflammaed and non-inflammed conditions.

Again I am using 0.1 as the adjusted p-value cutoff.

First let's examine the TCT model.

# Non-inflammed condition
tct_ni_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_NI_MKI vs NI_V_transcript_level.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

tct_ni_MK2i_DGE <- tct_ni_MK2i$X1[tct_ni_MK2i$padj < 0.1]

# Inflammed condition
tct_i_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_I_MKI vs I_V_transcript_level.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

tct_i_MK2i_DGE <- tct_i_MK2i$X1[tct_i_MK2i$padj < 0.1]

# Examine the overlap between NI and I conditions
tct_MK2i_overlap <- intersect(tct_i_MK2i_DGE, tct_ni_MK2i_DGE)

grid.newpage()
draw.pairwise.venn(length(tct_i_MK2i_DGE),
                   length(tct_ni_MK2i_DGE),
                   length(tct_MK2i_overlap)-1,
                   catergory <- c("MK2i target transcripts\nin inflammed TCT",
                                  "MK2i target transcripts\nin non-inflammed TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5750], polygon[GRID.polygon.5751], polygon[GRID.polygon.5752], polygon[GRID.polygon.5753], text[GRID.text.5754], text[GRID.text.5755], text[GRID.text.5756], lines[GRID.lines.5757], text[GRID.text.5758], text[GRID.text.5759])

Now let's examine the TNF model.

# Non-inflammed condition
tnf_ni_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_NI_MKI vs NI_V_transcript_level.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

tnf_ni_MK2i_DGE <- tnf_ni_MK2i$X1[tnf_ni_MK2i$padj < 0.1]

# Inflammed condition
tnf_i_MK2i <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_I_MKI vs I_V_transcript_level.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

tnf_i_MK2i_DGE <- tnf_i_MK2i$X1[tnf_i_MK2i$padj < 0.1]

# Examine the overlap between NI and I conditions
tnf_MK2i_overlap <- intersect(tnf_i_MK2i_DGE, tnf_ni_MK2i_DGE)

grid.newpage()
draw.pairwise.venn(length(tnf_i_MK2i_DGE),
                   length(tnf_ni_MK2i_DGE),
                   length(tnf_MK2i_overlap)-1,
                   catergory <- c("MK2i target transcripts\nin inflammed TNF",
                                  "MK2i target transcripts\nin non-inflammed TNF"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.05,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5760], polygon[GRID.polygon.5761], polygon[GRID.polygon.5762], polygon[GRID.polygon.5763], text[GRID.text.5764], text[GRID.text.5765], text[GRID.text.5766], lines[GRID.lines.5767], text[GRID.text.5768], text[GRID.text.5769])

Now let's take a looks and see if there is any overlap between TCT and TNF models.

# Inflammed conditions
i_MK2i_overlap <- intersect(tnf_i_MK2i_DGE, tct_i_MK2i_DGE)
grid.newpage()
draw.pairwise.venn(length(tnf_i_MK2i_DGE),
                   length(tct_i_MK2i_DGE),
                   length(i_MK2i_overlap)-1,
                   catergory <- c("MK2i target transcripts\nin inflammed TNF",
                                  "MK2i target transcripts\nin inflammed TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5770], polygon[GRID.polygon.5771], polygon[GRID.polygon.5772], polygon[GRID.polygon.5773], text[GRID.text.5774], text[GRID.text.5775], text[GRID.text.5776], lines[GRID.lines.5777], text[GRID.text.5778], text[GRID.text.5779])

# Non-inflammed conditions
ni_MK2i_overlap <- intersect(tnf_ni_MK2i_DGE, tct_ni_MK2i_DGE)
grid.newpage()
draw.pairwise.venn(length(tnf_ni_MK2i_DGE),
                   length(tct_ni_MK2i_DGE),
                   length(ni_MK2i_overlap)-1,
                   catergory <- c("MK2i target transcripts\nin non-inflammed TNF",
                                  "MK2i target transcripts\nin non-inflammed TCT"),
                   lty = "blank",
                   ex.text = FALSE,
                   fill = c("pink", "lightblue"),
                   cat.pos = c(200, 135), cat.dist = 0.08, margin = 0.1,
                   fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5780], polygon[GRID.polygon.5781], polygon[GRID.polygon.5782], polygon[GRID.polygon.5783], text[GRID.text.5784], text[GRID.text.5785], text[GRID.text.5786], lines[GRID.lines.5787], text[GRID.text.5788], text[GRID.text.5789])

Identify Biomarkers of effective response to MK2i.

Biomarkers are defined as genes that are targets of MK2i and are dysregulated in only TCT model of inflammation but not TNF model of inflammation.

First we will be more stringent, I will only use targets of MK2i in inflammed TCT model.

tct_biomarker_raw <- intersect(tct_i_MK2i_DGE, tct_inf_DGE)

tct_biomarker <- setdiff(tct_biomarker_raw, tnf_inf_DGE)

n23 <- intersect(tct_i_MK2i_DGE, tnf_inf_DGE)

n123 <- intersect(inf_overlap, tct_i_MK2i_DGE)

grid.newpage()
draw.triple.venn(area1 = length(tct_inf_DGE),
                 area2 = length(tnf_inf_DGE),
                 area3 = length(tct_i_MK2i_DGE),
                 n12 = length(inf_overlap)-1,
                 n23 = length(n23)-1,
                 n13 = length(tct_biomarker_raw)-1,
                 n123 = length(n123)-1,
                 catergory <- c("TCT inflammation DGE transcripts",
                                "TNF inflammation DGE transcripts",
                                "MK2i target transcripts\nin inflammed TCT"),
                 lty = "blank",
                 ex.text = FALSE,
                 fill = c("pink", "lightblue","lightgreen"),
                 cat.dist = 0.08, margin = 0.1,
                 fontfamily = "sans", cat.fontfamily = "sans")

## (polygon[GRID.polygon.5790], polygon[GRID.polygon.5791], polygon[GRID.polygon.5792], polygon[GRID.polygon.5793], polygon[GRID.polygon.5794], polygon[GRID.polygon.5795], text[GRID.text.5796], text[GRID.text.5797], text[GRID.text.5798], text[GRID.text.5799], text[GRID.text.5800], text[GRID.text.5801], text[GRID.text.5802], text[GRID.text.5803], text[GRID.text.5804], text[GRID.text.5805])

So there are three genes that popped out as potential biomarkers. List them here:

tct_biomarker

##  [1] "ENSMUST00000034586" "ENSMUST00000072646" "ENSMUST00000099960"
##  [4] "ENSMUST00000111224" "ENSMUST00000139385" "ENSMUST00000149654"
##  [7] "ENSMUST00000161561" "ENSMUST00000163854" "ENSMUST00000184716"
## [10] "ENSMUST00000187521"

We want to make sure the LFC are the opposite for these three genes in TCT inflammation DGE list and MK2i target in inflammed TCT.

tct_inf_DGE_LFC <- tct_inf$log2FoldChange[tct_inf$X1 %in% tct_biomarker]
tct_i_MK2i_LFC <- tct_i_MK2i$log2FoldChange[tct_i_MK2i$X1 %in% tct_biomarker]
biomarker_table <- cbind(tct_biomarker, tct_inf_DGE_LFC, tct_i_MK2i_LFC)
colnames(biomarker_table) <- c("Ensembl_ID", "LFC(TCT_inflammation_related)", "LFC(TCT_inflammed_MK2i_target)")
biomarker_table

##       Ensembl_ID           LFC(TCT_inflammation_related)
##  [1,] "ENSMUST00000034586" "1.55386067789831"           
##  [2,] "ENSMUST00000072646" "1.52312489295257"           
##  [3,] "ENSMUST00000099960" "2.99948049251717"           
##  [4,] "ENSMUST00000111224" "2.71392400685194"           
##  [5,] "ENSMUST00000139385" "1.99712975899567"           
##  [6,] "ENSMUST00000149654" "4.01314578122636"           
##  [7,] "ENSMUST00000161561" "-1.89168738340235"          
##  [8,] "ENSMUST00000163854" "5.26859392363033"           
##  [9,] "ENSMUST00000184716" "-1.4122366120438"           
## [10,] "ENSMUST00000187521" "-0.73802744116672"          
##       LFC(TCT_inflammed_MK2i_target)
##  [1,] "-1.176873529389"             
##  [2,] "-1.01910347165183"           
##  [3,] "-2.21805117040699"           
##  [4,] "-1.53677374422883"           
##  [5,] "-0.716209541424344"          
##  [6,] "-3.13359535848912"           
##  [7,] "-1.02593483398839"           
##  [8,] "-1.47752839171476"           
##  [9,] "-1.11455436327322"           
## [10,] "-0.430562088954025"

It appears that MK2i are downregulating 7/10 transcripts, while inflammation upregulates them in the TCT model. This indicates that they could be the genes that mediate the responsiveness to MK2i treatment.